library(librarian)
library(clusterProfiler)
simpl_enrich <- simplify
shelf(tidyverse,
      ggrepel,
      patchwork,
      ggpubr)
source('00_util_scripts/mod_bplot.R')

# examine mva pathway ----------
kegg_mva <-
  c('Hmgcr','Hmgcs1','Hmgcs2','Fdps','Mvd','Idi1','Idi2', 'Pmvk','Acat1','Acat2','Mvk') |>
  str_to_upper()

# examine interested cytokines ------
key_cytokine <- c('Il6','Ccl20','Tslp','Flt3lg','Csf2','Tnf') |>
  str_to_upper()

kera_slevhc <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_kera_sle-hc_deg.csv')

sle_v3hvl <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_sle_v3hvl.csv')

hc.v3hvl <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_hc_v3hvl.deg.csv')

v3h_slevhc <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_v3h_sle-hc_deg.csv')

v3l_slevhc <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_v3l_sle-hc_deg.csv')

### FIG: kera SLEvsHC volcano ============
kera_slevhc |>
  mutate(p_val_adj = ifelse(gene == 'CCL20', 1e-295, p_val_adj)) |>
  filter(p_val_adj != 0) |>
  plot_pub_volc('SLE', highlights = c('CCL20','TNF','FLT3LG','TRPV3'),force = TRUE)

publish_pdf('mission/FPP/figures/SLE-HC_kera_volcano.pdf', width = 65)

kera_slevhc_cyto <- kera_slevhc |>
  filter(gene %in% key_cytokine)

### FIG: up go kera SLE vs kera HC =======
kera_sle_upgo <- kera_slevhc |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T) |>
  simpl_enrich()

kera_sle_upgo@result |>
  head(10) |>
  publish_enrichment()

ggsave('figures/kera_SLE-HC_upgo.pdf',
       width = 65, height = 50, units = 'mm')

kera_slevhc |>
  right_join(chole_p_tib, join_by(gene == chole_p_str))

### FIG: HC v3h vs v3l volcano ==========
hc.v3hvl |>
  filter(p_val_adj != 0) |>
  plot_bill_volc('HC V3-hi KC', 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva)) +
  ggtitle('HC TRPV3-hi KC vs HC TRPV3-lo KC')

hc.v3hvl |>
  filter(p_val_adj != 0) |>
  plot_pub_volc('HC V3-hi KC', 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/HC_v3h-v3l_volcano.pdf')

### FIG: HC v3h vs v3l upgo ==========
hc.v3hvl_upgo <- hc.v3hvl |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

hc.v3hvl_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in HC V3-hi KC vs HC V3-lo KC')

hc.v3hvl_upgo@result |>
  head(10) |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in HC V3-hi KC vs HC V3-lo KC')

publish_pdf('mission/FPP/figures/HC_v3h-v3l_upgo.pdf', width = 65)

hc.v3hvl_upgo <- hc.v3hvl |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'ALL',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

hc.v3hvl_upgo@result |>
  write_csv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.upgo.csv')

read_tsv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.go.digest.tsv') |>
  plot_enrichment(n = 12) +
  labs(title = 'Upregulated GO pathways in human HC v3h vs v3l KC')

### HC v3h vs v3l kegg =========
hc.v3hvl_upkegg <- hc.v3hvl |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(gene) |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
       OrgDb = 'org.Hs.eg.db') |>
  pull(ENTREZID) |>
  enrichKEGG(qvalueCutoff = .1)

hc.v3hvl_upkegg@result |>
  filter(qvalue < .1) |>
  write_csv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.upkegg.csv')

hc.v3hvl_upkegg@result |>
  filter(qvalue < .1) |>
  plot_enrichment() +
  labs(title = 'Upregulated KEGG pathways in human HC v3h vs v3l KC')

### HC v3h vs v3l: GSEA ===========
hc.v3hvl_gse.kegg <- hc.v3hvl |>
  pull(gene) |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
       OrgDb = 'org.Hs.eg.db') |>
  left_join(hc.v3hvl, join_by(SYMBOL == gene)) |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T) |>
  gseKEGG()

hc.v3hvl_gse.kegg@result |>
  DT::datatable()

### FIG: SLE v3h vs v3l volcano =========
sle_v3hvl |>
  filter(p_val_adj != 0) |>
  plot_bill_volc('SLE V3-hi KC', 'SLE V3-lo KC',
                highlights = c(key_cytokine, kegg_mva))

sle_v3hvl |>
  filter(p_val_adj != 0) |>
  plot_pub_volc('SLE V3-hi KC', 'SLE V3-lo KC',
                highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/SLE_v3h-v3l_volc.pdf')

### FIG: SLE v3h vs v3l upgo ==========
sle_v3hvl_upgo <- sle_v3hvl |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

sle_v3hvl_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in SLE V3-hi KC vs V3-lo KC')
  
sle_v3hvl_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in SLE: V3-hi KC vs V3-lo KC')

publish_pdf('mission/FPP/figures/SLE_v3h-v3l_upgo.pdf', width = 65)

## FIG: v3h SLE-HC volcano ===========
v3h_slevhc |>
  filter(p_val_adj > 0) |>
  plot_bill_volc(group1 = 'SLE V3-hi KC', group2 = 'HC V3-hi KC',
                highlights = c(key_cytokine, kegg_mva))

v3h_slevhc |>
  filter(p_val_adj > 0) |>
  plot_pub_volc(group1 = 'SLE V3-hi KC', group2 = 'HC V3-hi KC',
                highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/v3h_SLE-HC_volcano.pdf')

## FIG: v3h SLE-HC upgo ========
v3h_slevhc_upgo <- v3h_slevhc |>
  filter(p_val_adj < .05 & avg_log2FC > 1) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

v3h_slevhc_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in SLE V3-hi KC vs HC V3-hi KC')

v3h_slevhc_upgo@result |>
  head(10) |>
  publish_enrichment()

publish_pdf('figures/v3h_SLE-HC_upgo.pdf', width = 60)

## FIG: v3lo SLE-HC volcano ==========
v3l_slevhc |>
  filter(p_val_adj > 0) |>
  plot_bill_volc(group1 = 'SLE V3-lo KC', group2 = 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

v3l_slevhc |>
  filter(p_val_adj > 0) |>
  plot_pub_volc(group1 = 'SLE V3-lo KC', group2 = 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/v3l_SLEvHC_volcano.pdf')

## FIG: v3lo SLE-HC upgo ========
v3l_slevhc_upgo <- v3l_slevhc |>
  filter(p_val_adj < .05 & avg_log2FC > 1) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

v3l_slevhc_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in SLE V3-lo KC vs HC V3-lo KC')

v3l_slevhc_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in V3-lo KC: SLE vs HC')

publish_pdf('mission/FPP/figures/v3l_SLE-HC_upgo.pdf', width = 65)
